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We present a novel method, which we shall refer to as the dual minima hopping method (DMHM), 
that allows us to find the global minimum of the potential energy surface (PES) within density 
functional theory (DFT) for systems where a fast but less accurate calculation of the PES is possible. 
This method can rapidly find the ground state configuration of clusters and other complex systems 
with present day computer power by performing a systematic search. We apply the new method 
to silicon clusters. Even though these systems have already been extensively studied by other 
methods, we find new configurations for Sii6, Siis and Siig that are lower in energy than those 
found previously. 



Determining the structure of a molecule, cluster or 
crystal is one of the most fundamental and important 
tasks in solid state physics and chemistry. Practically all 
physical properties of a system depend on its structure. 
The structural configurations of a system are determined 
by the Born-Oppenheimer PES, which gives the energy of 
a system as a function of its atomic coordinates. Minima 
of the PES give stable configurations. The global mini- 
mum gives the ground state configuration. At low enough 
temperature the system will be found in this global mini- 
mum structure assuming that this structure is kinetically 
accessible. Since the zero point energy of different struc- 
tures varies negligibly, the determination of the ground 
state structure is equivalent to the mathematical problem 
of finding the global minimum of the PES. 

It is well established that the PES of a condensed mat- 
ter system can be calculated with good accuracy within 
DFT. Nevertheless, DFT methods have not been used 
up to now as a standard tool in algorithms that attempt 
to determine the ground state of complex systems be- 
cause most algorithms for the determination of the global 
minimum require an enormous number of evaluations of 
the PES. Since each evaluation requires a full electronic 
structure calculation, these algorithms are computation- 
ally too demanding within the full DFT framework. A 
systematic search for the global minimum is however pos- 
sible with cheaper methods such as tight binding and 
force field methods. 

In summary, with present methods one has either the 
choice of using methods with a limited power of pre- 
dictability or of doing a constrained search for the global 
minimum. In a constrained search one fixes some atomic 
positions or imposes some structural motifs, but experi- 
ence shows that the global minimum is often missed in 



this way. To overcome this dilemma several researchers 
have adopted an approach where one first effectuates a 
systematic search with a method that allows for a fast 
but inaccurate calculation of the PES to obtain some 
candidate structures. Which of the candidate structures 
is lowest in energy is determined in a second step by 
DFT calculations. As we shall show later this approach 
is generally not applicable. 

Other researches have combined systematic search al- 
gorithms with DFT methods, but their algorithms re- 
quired too many DFT calculations to be computation- 
ally feasible if one wants to find the global minimum. 
Rothlisberger et al Q have used simulated annealing 
within DFT to find structural motifs of the mid-size clus- 
ters but their final lowest energy geometries were ob- 
tained by other means. Yoo and Zeng 0, have 
combined basin-hopping (BH) with DFT and were able 
to find new low- lying minima for some clusters, among 
them Sii6, Sii7 and Siis- For Sii6, they have found a 
new global minimum structure by performing a system- 
atic BH search within DFT. As we shall see, both the 
systematic BH for Sii6 as well as the constrained BH for 
Sii7 and Siis within DFT have missed the global mini- 
mum. 

In this paper we shall present a method that allows for 
a systematic search for the global minimum of the PES of 
a complex system within DFT. The method is a modifi- 
cation of the minima hopping method (MHM) |4( . In the 
MHM one visits a series of local minima until the global 
minimum is found. The algorithm has a double loop 
structure. In the inner loop one attempts to escape from 
the current minimum, in the outer loop one accepts or 
rejects new minima found by successful escape attempts. 
A history list keeps track of all minima found. A feed- 
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back mechanism uses information from this history list to 
make more vigorous escape attempts when the algorithm 
is revisiting previously found minima thereby preventing 
the algorithm from getting trapped in an incorrect min- 
imum. The inner escape loop contains two basic steps. 
The first does a certain number of molecular dynamics 
(MD) moves until one has overcome at least one energy 
barrier. The second step consists in performing a stan- 
dard geometry relaxation to reach the closest minimum 
with an accurate method. 

In the ordinary version of the MHM the forces for 
the MD and for the geometry optimization part are done 
with the same method. Fast methods such as force field 
or tight binding methods have to be used to limit the 
computing time to an acceptable length. In the modified 
MHM presented in this paper two different methods are 
combined: a slow but accurate method and a fast but 
less accurate method. The fast method is used for the 
MD part and for the first few steps of the geometry opti- 
mization. The accurate method is then used for the final 
geometry optimization and the evaluation of the energy 
of the relaxed structure. In this way the search for the 
global minimum is reduced to a relatively small number 
of geometry optimizations with the accurate and expen- 
sive method plus a much larger number of force evalua- 
tions with the fast method. Henceforth, we shall refer to 
this modified minima hopping algorithm, that combines 
the two methods for the calculation of the forces, as the 
dual minima hopping method (DMHM). 

The fact that the input configuration for the geome- 
try optimization with the accurate method is a configu- 
ration that was prerelaxed with the fast method is im- 
portant for the stability of the entire algorithm if the 
accurate method is a DFT method. DFT programs do 
typically not converge if the input configuration is far 
from any physically reasonable configuration. The pre- 
relaxation with the fast method excludes the possibility 
that a physically unreasonable state is used as an input 
configuration. From the previous considerations it might 
seem advantageous to do a full prerelaxation, i.e. to use 
a minimum of the fast method as the input for the geom- 
etry optimization with the accurate method. If the fast 
method is a reasonable approximation then a local min- 
imum found by it will often be close to a local minimum 
of the accurate method. 

Unfortunately, in general there is no one-to-one corre- 
spondence between minima obtained from the two meth- 
ods. Therefore, some minima obtained using the accu- 
rate method are inaccessible from the starting configura- 
tions provided by the fast method. For this reason only 
a small number of steps should be done in the prerelax- 
ation with the fast method. In this way the ensemble of 
the starting configurations for the geometry optimization 
with the accurate method comprises a considerable part 
of the configurational space (and not only the ensemble 
of all the minima of the fast method) and one can reach 



virtually any minimum of the accurate method. 
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FIG. 1: The truncated polynomial approximation of the 
Lennard-Jones potential and the exact Lennard-Jones poten- 
tial. 

The Bell-Evans-Polanyi (BEP) principle Q states that 
highly exothermic chemical reactions have a low activa- 
tion energy. In the context of a global minimum search, 
this means energetically low configurations will be prefer- 
ably found behind low energy barriers. The BEP princi- 
ple is essential for the success of the MHM as has been 
shown in 4] . The correlation between the barrier height 
and the energy of the minimum 'behind' the barrier cer- 
tainly deteriorates if one is combining two different meth- 
ods. This implies more local minima will be visited, on 
average, with the DMHM before the global minimum is 
found than with the ordinary MHM. In order to explore 
the influence of this reduced correlation we did system- 
atic tests with a 38 atom Lennard-Jones (LJ) cluster. 
This is a system for which the global minimum is hard to 
find since it is contained in a small secondary funnel 0, 
but the computing time is small since the potential can 
be evaluated very rapidly. As the accurate method we 
used the LJ potential. As the 'fast' method we used a 
truncated polynomial approximation of the LJ potential 
as shown in Fig.Q As expected, the number of local min- 
ima that are visited on average before the global mini- 
mum is found increases from 380 to 530, nevertheless, the 
number of force evaluations needed with the 'expensive' 
exact LJ method is reduced by a factor of 5. 

To demonstrate that the DMHM can indeed find the 
global ground state geometry of real clusters, we have ap- 
plied it to silicon clusters. Numerous groups are involved 
in the search of the ground state of silicon clusters and 
there are at least 50 theoretical papers on this subject Q- 
0- Applying DMHM to silicon clusters we were 
able to find within several days of computing time all of 
the known structures Q , ^} , , m the range Si4- 
Siig and we even found lower energy structures for Sii6, 
Sii8 and Siig in spite of the fact that silicon clusters up to 
19 atoms in size have already been extensively studied. 
The new global minimum structures within CPMD/PBE 
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(see below) Sii6«, Siisa and Siig a as well as the new low- 
lying isomers Sii6&, Sii7 a and Siub are shown in Fig. El 
The structure Sii6a contains the TTP-Sig-subunit |l8q 
and is compact in contrary to the structure Sii6 reported 
by X.C.Zeng |3|. The structure Sii8 a is prolate and con- 
sists of two TTP-Sig-subunits which are rotated against 
each other. The structure Siig a consists of a TTP-Sig- 
subunit and a Siio-subunit. The low- lying isomer Sii66 
is compact and highly symmetric. The low- lying isomer 
Sii7 a consists of a TTP-Sig-subunit and a Sis-subunit. 
The low-lying isomer Siub consists of two equal 7-blocks, 
which are rotated against each other, and a triangle as 
a cleaving block. In contrast to the previous works, our 
configurations were found by the DMHM automatically 
after having visited only a few hundred local minima. 
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FIG. 2: Lowest energy geometries Sii6a, Sii66, Sii7a, Sim, 
Siisa and Siig a found in this work with DMHM and the pu- 
tative global minimum structures Sii6 3], Siis E3 an d 
Siig [13 reported previously. The new geometries will be 
posted on the Cambridge Cluster Database [T^. 



As the fast method we have used the Lenosky tight 
binding scheme for silicon [2£|. The accurate method is 
DFT as implemented in the Quickstep code 21]. Af- 
ter having performed the DMHM with Quickstep using 
a relatively small Gaussian basis set and the local den- 
sity approximation (LDA), we have calculated accurate 
final energies and zero-point energies with the CPMD 
program [23 using the PBE functional j23| , a high accu- 
racy pseudo-potential Q, large super-cells (24 A) and 
a sufficient plane wave cutoff (28 Rydberg). The results 
for Sii6a, Siisa and Siig a as compared to Sii6, Siis and 
Siig are presented in Table |I] The low- lying isomer Sii66 
is virtually isoenergetic with the structure Sii6- The iso- 
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-0.15 


-0.01 


-0.08 


PBE(Z) 


-0.16 


-0.07 


-0.09 



TABLE I: The energy differences in eV without and with 
zero-point energy correction between the lowest energy ge- 
ometries Sii6a, Siisa and Siig a found in this work with 
DMHM and the putative global minimum structures Sii6 0, 
oils and Siig reported previously using the PBE 

exchange- correlation functional as implemented in CPMD. 



mer Sii7 a is lower by 0.16 eV, the isomer Si^ is lower by 
0.06 eV within CPMD/PBE than the new low-lying Sii 7 
isomer reported by X.C.Zeng 3]. The Sii 7 structure re- 
ported by Ho et al. in 7] is however 0.10 eV lower than 
our isomer Sii 7a . In contrast to other exchange correla- 
tion functionals, the PBE functional [2^ was not fitted 
to any chemical systems with simple bond structures and 
is expected to give the most accurate description of the 
complex bonding patterns found in silicon clusters. The 
term 'accurate' must be handled with caution however, 
since DFT is only an approximation and, as a matter of 
fact, the energetic ordering may change if one uses dif- 
ferent functionals [l^ . 

Among the various force fields and tight binding 
schemes we have tested, the Lenosky tight binding 
scheme |2£j gave the best agreement with the DFT en- 
ergies. It can predict the DFT energies with an error of 
roughly 1 eV as shown in Fig. 03 Fig. also shows why 
the common approach of first finding candidate struc- 
tures by doing a systematic search with a cheap method 
and then checking by an accurate method which of the 
candidate structures gives the global minimum is prob- 
lematic except for very small systems. For a 25 atom 
silicon cluster the number of geometric configurations 
within 1 eV above the ground state is of the order of 
10 4 states, for a 33 atom cluster it is already of the order 
of 10 5 states and it increases exponentially with system 
size. It is therefore virtually impossible to check which 
out of these 10 4 to 10 5 configurations is the global min- 
imum in DFT. Besides, because of the absence of the 
one-to-one correspondence between the local minima of 
the fast method and of the accurate method, it is not 
guaranteed that any of the minima of the fast method 
will lead to the global minimum of the accurate method 
upon relaxation. 

The identification of the previously visited minima is 
an essential ingredient of the MHM. In the context of 
the ordinary MHM the energy can be used to identify 
configurations since it is possible to calculate the energy 
with many significant digits both for force fields and tight 
binding schemes. With DFT programs this is not any 
more possible because of the presence of numerical noise. 
For this reason we have used in addition to the energy all 
inter- atomic distances. Two DFT minima are considered 
to be identical if all their inter-atomic distances ordered 
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FIG. 3: The correlation between tight binding and density 
functional energies for various configurations of a Si25 cluster. 
If the correlation was perfect all the points would lie on the 
diagonal. Instead the scattering shows that the tight binding 
energies can predict the energy differences between various 
cluster configurations only with an error of about 1 eV. 



by magnitude agree to within a certain tolerance. 

In summary, we have presented a method that allows 
one to find the global minimum of the DFT potential en- 
ergy surface within acceptable computer time for mod- 
erately complex systems. The method is efficient for the 
following reasons. First, it requires only DFT calcula- 
tions for configurations where DFT programs typically 
converge without problems. It does not, for instance, 
require DFT calculations for configurations generated 
by random displacements from a previous configuration. 
Second, the MHM is highly efficient in the sense that 
the number of minima visited before the ground state is 
found is small. Even though the DMHM is not quite as 
good from this perspective it is still efficient if the fast 
method used for the MD part is qualitatively correct. 
Third, most of the force evaluations are done with the 
fast method and the total effort for finding the global 
minimum is equal to the effort of doing only an afford- 
able number of geometry optimizations with the accurate 
method. 
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